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Abstract 

Interfacial mass transfer of low-diffusive substances in an unsteady flow environment is marked by a very thin 
boundary layer at the interface and other regions with steep concentration gradients. A numerical scheme capable 
of resolving accurately most details of this process is presented. In this scheme, the fourth-order accurate WENO 
method developed by Liu et al. (1994) was implemented on a non-uniform staggered mesh to discretize the scalar 
convection while for the scalar diffusion a fourth-order accurate central discretization was employed. The discretiza- 
tion of the scalar convection-diffusion equation was combined with a fourth-order Navier-Stokes solver which solves 
the incompressible flow. A dual meshing strategy was employed, in which the scalar was solved on a finer mesh than 
the incompressible flow. The solver was tested by performing a number of two-dimensional simulations of an unstably 
stratified flow with low diffusivity scalar transport. The unstable stratification led to buoyant convection which was 
modelled using a Boussinesq approximation with a linear relationship between flow temperature and density. The or- 
der of accuracy for one-dimensional scalar transport on a stretched and uniform grid was also tested. The results show 
that for the method presented above a relatively coarse mesh is sufficient to accurately describe the fluid flow, while 
the use of a refined mesh for the low-diffusive scalars is found to be beneficial in order to obtain a highly accurate 
resolution with negligible numerical diffusion. 
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1. Introduction 

To accurately resolve low-diffusivity scalar transport problems special numerical schemes are necessary for the 
discretization of the convective term in order to avoid under- and/or overshoots of the scalar quantity. The 
first order upwind method, for example, could be used to effectively avoid such under- and/or overshoots but 
at the cost of introducing an excessive amount of numerical diffusion dill- Up to now, a number of DNS studies of 
gas transfer across the air- water interfaces have been carried out for shear driven and stirred vessels. Hasegawa and 
Kasagi (3]|4) studied wind- shear driven mass transfer across the turbulent interface at a Schmidt number of Sc = 100. 
They used a pseudo-spectral Fourier method for the spatial discretization in the horizontal directions, whereas the 
finite volume method is employed in the normal direction in which turbulent and molecular mass fluxes are evaluated 
at a cell surface with second-order accuracy. Handler et al. [5 ] used a pseudo spectral approach with Fourier expan- 
sions to carry out direct numerical simulations for the transport of a passive scalar at a shear-free boundary in fully 
developed channel flow. Similarly, Banerjee also used a pseudo- spectral method to extensively study the mechanisms 
of turbulence and scalar exchange at the air- water interface in several publications (see (61I3 and references therein). 
Schwertfirm and Manhardt (8) also studied passive scalar transport in an turbulent channel flow for Schmidt numbers 
up to Sc = 49. They used a similar approach to the present work by solving the scalar on a finer grid then the ve- 
locity which was mapped by a conservative interpolation to the fine-grid. An explicit iterative finite- volume scheme 
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of sixth-order accuracy was employed to calculate all convective and diffusive fluxes, while for the time-integration a 
third order Runge-Kutta method was used (9). 

The pseudo-spectral methods used above have excellent error properties when the solution is relatively smooth. 
However, a main disadvantage of using spectral methods lies in the formation of non-physical oscillations near steep 
gradients that may regularly occur in the solution of a convection-diffusion problem when the diffusivity is extremely 
small. These Gibbs oscillations (under- and/or overshoots) ifTOli near steep gradients are not uncommon and can 
also be found when higher-order central finite-difference methods are used on a relatively coarse mesh to discretize 
the convection of a scalar with low diffusivity. This can be overcome by using weighted essentially non-oscillatory 
schemes (WENO) which have excellent shock capturing capabilities. Their non-oscillatory behaviour is advantageous 
in dealing with very steep gradients. In this paper, we present a numerical scheme developed specifically to resolve the 
details of the interfacial mass transfer of low-diffusive substances by adapting the weighted essentially non-oscillatory 
(WENO) by Liu et al. ifTTTl . They are based on essentially non-oscillatory (ENO) schemes which were first published 
in the meanwhile classic paper of Harten et al. fl2l . Liu et al. IfTTTl introduced the idea of taking a convex combination 
of interpolation polynomials to construct a stencil using non-linear weights with a high order-of-accuracy in smooth 
regions while weighing out the non- smooth stencils in regions containing steep gradients or discontinuities. They 
studied WENO(2r - 1) schemes for different stencil sizes, i.e. r = 2 (WEN03) and r = 3 (WEN05). 

In the meantime a large variety of WENO schemes has been developed. Many improvements were made by 
modifying the smoothness determination. For instance, Jiang and Shu fl"3l introduced a new smoothness indicator 
that is used to evaluate the non-linear weights. The size of the stencil has also been further increased by Balsara and 
Shu lfT4l extending it up to r = 6 (WENOll). Henrick et al. [15 ] could show that the weights generated by the 
classical choice of smoothness indicators in Hl3l failed to recover the maximum order of the scheme at critical points 
of the solution where the first derivatives are zero. They developed the so called WENOM schemes where a mapping 
procedure is introduced to keep the weights of the stencils as close as possible to the optimal weights. The resulting 
(mapped) WENOM scheme of 031 presented more accurate results close to discontinuities. Even more recently, 
Borges et al. [16] achieved the same results as mapped WENO schemes without mapping but by improving the 
accuracy of the classical WEN05 scheme by devising a new smoothness indicator and non-linear weights using the 
whole 5-points stencil and not the classical smoothness indicator of Jiang and Shu [13 ] which uses a composition of 
three 3-points stencils. The schemes of Borges et al. are known as WENO-Z schemes. Of all the schemes discussed 
above the classical WEN05 scheme is used most widely (15J [T71 [THJ. In our simulations, we do not expect any 
discontinuities in the scalar field so that the classical WEN05 scheme of Liu et al. ifTTIl is a good choice to accurately 
resolve low-diffusive scalar transport which may lead to steep concentration gradients. 

An example of application is given for the 2D case of buoyant-convectively driven mass transfer with a Prandtl 
number of Pr = 6 and Schmidt numbers of S c = 20, 50, 200 and 500. One typical process in nature of such a case 
is the absorption of oxygen into lakes during night time. This process is controlled by the low diffusivity of the 
dissolved gas in the water and the convective instability triggered by the density difference between the cold water at 
the top surface and the warm water in the bulk. The convective-instability enhances the gas transfer into the water 
body significantly compared to the static condition with only diffusive gas transfer. This low diffusive process results 
in a very thin concentration boundary layer that is found at the water surface. Experimental measurements near the 
surface (such as the mass flux) are very difficult. There is a need to fully resolve the near surface mechanisms in 
order to understand the physical mechanism. Below, the capability of the newly developed code to accurately resolve 
low-diffusivity scalar transport problems will be illustrated. 

2. Formulation of Numerical Method 

2.7. Discretization of the convection-diffusion equation of the scalar (p 

The scheme that is employed here is a variant of the WEN05 finite difference scheme as described by Liu et al. 
AH . The WENO schemes use an approximation of the scalar fluxes at the cell interface by employing interpola- 
tion schemes. The reconstruction procedure produces a high order accurate approximation to the solution from the 
calculated cell averages. 

In this section we briefly outline the discretization of the transport equation for the scalar (p. The two-dimensional 
convection diffusion equation of (p = (p(x, z, t) reads 
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where x and z are the horizontal and vertical directions, respectively and t denotes time. The diffusive term on the right 
is discretized using a fourth-order accurate central scheme, while the convective term is discretized using a variant of 
the fourth-order WENO scheme developed by Liu et al. |[TT1 . Below the implemented scheme is detailed only in one 
dimension. Generalization to higher dimensions is straightforward. 

When ignoring the diffusive term, the one dimensional variant of ([T]) can be rewritten as 

— = -u— (2) 
dt dx' 

where u is the velocity in the x-direction. As we employ a staggered mesh, for the volume centred around x = Xi, the 
convective fluxes Rt and Rt are defined by 

Rt = Pi(* i+ h)+ PM(*i + 0) 

^0 + ^1+^2 2 ®0 + ®l + ®2 2 + d\ + Cl2 2 

and 

R- = Pi-i(x i i)+ ^ Pi + i(Xi-i) (4) 

<2q + a\ + Cl2 2 <2q + a\ + Cl2 2 <2q + <2i + <22 2 

where «o>^i>^2 are coefficients and P the Lagrange interpolations polynomials defined in ([6]). In the above, the 
weights for the convex combination of the quadratic Lagrange interpolation polynomials are given by, 
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where s = 10 9 and the smoothness indicator ISi is defined by 



ISi = ^(((fi-i - (fi-if + ((ft - (fii-i) 2 ) + (<Pi - 2(fi-i + (fi-if- 
The modified quadratic Lagrange interpolations Pi(x) in equations ([3} and ^ read 
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High-order polynomial interpolations to the midpoints x i+ \ are computed using known grid values of the scalar (p. 
The scheme uses a 5-points stencil which is divided into three 3-points stencils as shown in Fig.[T] 

These are interpolations of the scalar to the faces of the volume combined with a smoothing term at the right. 
Using the above, depending on the signs of u t _ i and u i+ 1 , we have four possible ways to calculate the discretization 

of the convective terms Li(<p) = (-m^) U m x f- 
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Figure 1: Schematic illustration of the weighted 5 point convex combination composed of three 3 -points stencils S ,Si,S 2 and 
their respective weights a ,ai,a 2 used in the classical WEN 05 scheme. l[T9\l . 
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The diffusive term on the right hand side of ([]} is discretized using a fourth-order central finite difference method 
for the second derivative such as, 

d 2 (f _ -<fi+2,k + ~ 30^ + ~ (fj-u m 
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and 
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where = - x t _i) and = - z t _i) , respectively. On a stretched mesh the actual discretisa- 

tion coefficients are obtained from the above equations using Lagrange interpolations. The time integration of the 
convection-diffusion equation is implemented using a third order Runge-Kutta method (RK3) developed by Shu and 
Osher lfT9l that reads, 



= ^ + |^ + |Atf,(^) (9) 



In the case of using three quadratic interpolations a fourth-order accuracy can be achieved. Note that the weights 
given to the interpolating polynomials depend on the local smoothness of the solution. Interpolation polynomials 
defined in regions where the solution is smooth are given higher weights than those in regions near discontinuities 
(shocks) or steep gradients (like the gas concentration near the interface in this case). 
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2.2. Flow Solver 

This section outlines the numerical method of the flow solver used in the buoyancy driven flow. The velocity 
field is solved by a finite-difference discretization of the convective terms using a fourth-order unconditionally kinetic 
energy conserving method combined with a fourth-order accurate central method for the diffusive terms l20l . The 2D 
incompressible Navier-Stokes equation is discretized on a non-uniform, staggered mesh in combination with a second- 
order accurate Adams-Bashforth time integration. The Boussinesq approximation is applied in order to account for 
the effects of buoyancy. For two-dimensional incompressible flow the continuity equation reads, 



du dw 
dx dz 



0, 



which in discretized form on a mesh as shown in Fig. 2a reads 
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The momentum equations are given by 
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where p is pressure and a and c represent the sum of the convective and diffusive terms 
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where Re = 100, based on a characteristic length scale of L = 1cm and a characteristic velocity ofu = w= lcm/ s The 
buoyancy term /3(T*) in equation fl5] ) is a function of the non-dimensional temperature T* (see ( [26] )) and is modeled 
using the Boussinesq approximation. In discretized form J3(T*) reads, 



0(T*) = -0.5611516 



1 k^ 1 k+l 



(16) 



The temperature difference between the top and the initial bulk temperature is only 3 °C. The temperature range 
in this particular case is between 20 and 23 °C where the relation between density and temperature can be assumed to 
be linear. 

When substituting the momentum equation into the continuity equation a Poisson equation for the pressure is 
obtained. The Poisson equation is iteratively solved using the conjugate gradient method with a diagonal precon- 
ditioning. From the obtained pressure field the new velocity field can be calculated by rearranging the discretized 
equations of fl2| ) and ( [T3] ), 



I n+l _ n+l \ 

<\ =w" , +Atl- P > Ml ~ P >* 1 +c" ,1. (18) 

Fig.|2a|shows the location of variables on a staggered mesh. To achieve kinetic energy conservation interpolations 
are required to evaluate the convective term. For instance at (x j+ i_,Zk) only the u- velocity component is available at 



u. 



i-l/2,k 



i,k+l/2 



O 



(T* p,cp) . . 



W i,k-l/2 
(a) Variables on a staggered mesh 



U 



i+l/2,k 



(T*ip) 



€>-- 

(T*p,cp) 



- O 

lT",<p) 



- o - 

IT'. 43) 



-^■O — O — r*-0 - 

(T*,cp) ; (T*,cp) ; fl*(p) 



(T*,ip) | (T* p,cp) | (T*,cp) 



■f-j t I I' 

-♦O —r*- O — f»>0 - 

(T*<p) : (T*,cp) i (T*,cp) 



(b) Dual sub-mesh refined by factor R = 2 (c) Dual sub-mesh refined by factor R = 3 



Figure 2: Variables on the new dual mesh. The flow field is solved on the outer coarse mesh, whilst the scalar is computed on a 
refined subgrid. For the transport of the scalar the velocities are interpolated onto the midpoints of the sub grid 

that location, while the w -velocity is only available at fez^+i). Hence an interpolation of w to the position where 
u is defined gives w i+ ^ k . An equivalent procedure for the z -momentum where u needs to be interpolated where w is 
defined gives u i k+ 1 . This yields to the discretization of the convective terms by a fourth-order central discretization, 



C x (u,w) i+ i k = - \ 
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The convective terms in the z -direction are discretized in a similar manner. The diffusive terms are discretized using 
the fourth-order accurate central discretization scheme (([7]) and ([8])) in which the coefficients on a non-uniform mesh 
in the z-direction are determined using Lagrange interpolations. 

2.3. Dual Mesh Approach 

Because the diffusivity of the scalars of interest is up to three orders smaller than that of the momentum, the 
resolution requirements for the flow field is less stringent as shown by the mesh refinement tests in Section]?] Thus, 
to save computing time a dual mesh approach is used as illustrated in Fig. [2] The velocity is solved on a coarser base 
mesh (Fig. [2a]), while the scalar is defined on the finer subgrid (Fig. [2b] and [2c]) so that the required computational 
resources are significantly reduced. 

To calculate the convective transport of the scalar the velocities are interpolated onto a subgrid using a fourth-order 
interpolation. When employing a subgrid refinement by a factor of R = 2 (Fig. [2b]) an interpolation is required for 
each subcell as the velocity location and the scalar locations do not coincide with their counter parts on the base mesh. 



In contrast Fig. [2c] shows that in the case of a subgrid refinement by a factor R = 3 some velocities and the central 
subcells for a scalar are defined at the same locations. 



3. ID Numerical Experiments 

In this section we apply the WENO-scheme for different test problems with the purpose to predict the accuracy 
of the method on uniform and stretched meshes, respectively. As the problem is a convection-diffusion problem 
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the numerical scheme was tested for both, scalar transport by convection only and by diffusion only. By using the 
quadratic Lagrange interpolation for reconstruction we expect to achieve fourth-order accuracy for the convective 
scalar transport. 

In both of the following cases, we use the previously described WENO scheme for spatial discretization and the 
3 rflf -order Runge-Kutta- scheme for time integration of the one-dimensional convection equation. Because we want to 
predict the numerical error in the WENO scheme, physical diffusion will be neglected. If (p(xu t) and if exact are the 
numerical and the exact solutions, respectively at (x t , t), the Li discretization error is given by: 

1 N 

^1 = — ^ I (<f(Xi, - (fexact) I 'SXi (20) 
i=l 

where N describes the number of nodes in the domain, t the time, i the node number and Sxt = \ (x i+ 1 - x t _ 1 ) 

3.1. Scalar transport by convection on uniform grid 

At first the WENO scheme is tested on uniform meshes. Therefore the following differential equation has been 
used to describe scalar convection within a one-dimensional domain: 



^ + *tt = 
dt dx 
0< x<2 



(21) 



using (p(x, 0) = (fo(x) and employing periodicity in space with period 2 for all t. 

The scalar distribution is initialized by a sine wave function (fio(x) = sin(7rx). Table [T] gives the Li error after 
running the simulation during one time-unit as well as the resulting order of accuracy. For this calculation a CFL- 
number of CFL = 0.5 has been used. 



N 


Li -error 


order 


40 


1.20 E-03 




80 


9.20 E-05 


3.7 


160 


4.30 E-06 


4.4 


320 


7.90 E-08 


5.8 



Table 1: Absolute error and order of convergence on uniform meshes. 



Starting from N = 40 nodes the Li error is decreasing when increasing the number of nodes to 80, 160 and 320. 
From N = 160 to N = 320 it can be seen that almost a 6 th order of accuracy is achieved. Although there are slight 
differences to Liu et al. [1 1 j, the implemented WENO-scheme reaches a high order of accuracy on uniform meshes. 

3.2. Scalar transport by convection on non-uniform meshes 

Using the modified Lagrange interpolations the WENO-scheme has been applied on non-uniform meshes on which 
the node distribution is given by: 



jc(0 = tan ((/ - 0.5 • N)h). < i < N (22) 
where h = ^ . The mesh is then scaled so that x(0) = and x(N) = 2. 
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N 


Li -error 


Li -order 


40 


3.70 E-03 




80 


5.30 E-04 


2.8 


160 


7.60 E-05 


2.8 


320 


1.04 E-05 


2.9 



Table 2: Absolute error and order of convergence on non-uniform meshes. 



The same conditions for calculating the errors and distribute the scalar in the domain have been used as for the 
uniform mesh. Due to the variable node distances Axi caused by the stretching, the CFL-number changes within the 
domain. However the maximal CFL-number has been set to CFL max = 0.5. 

Table [2] gives the results of the testing. Compared to uniform meshes it can be seen that the order of accuracy 
decreases to approximately 3. 

3.3. Scalar transport: pure diffusion 

In this section we applied the fourth-order accurate central discretization ([7]) for the solution of scalar diffusion. In 
the one-dimensional case, concentration gradients in the y- and z -direction are zero, and we have the one-dimensional 
diffusion equation for a scalar (p(x, t) 



i +D J = ° (23) 

where in the case of diffusive gas transfer into a liquid D = For the transfer of oxygen into water we 

have Sc = 500 and Re = 100, based on a characteristic length scale of L = 1cm and a characteristic velocity of 
u = lem/s. Hence the time unit is non-dimensionalised using L/U. A one dimensional domain has been chosen with 
<= x <= 5L. At x - the boundary condition ^(0, = 1 was imposed. A mesh with grid points n x was defined 
where the mesh was finer near the surface where the concentration boundary layer will form. The analytical solution 
for this boundary value problem is given by 



( p( x ,t) = l-erf\^=\ (24) 



The scalars (p are non-dimensionalized using the scalar magnitude <po at the top boundary and the magnitude in 
the bulk (f B so that 



(p = (25) 

<P0 ~ (fB 



while the temperature is made non-dimensional with, 



= T T 
Tb-T 

Figure [3] shows a comparison of numerical and theoretical results of the concentration boundary layer after t = 
80 L/U. Note that this plot only shows every 4 th gridpoint in the region near to the boundary where the concentration 
gradient occurs. 

The numerical results show very good agreement with the analytical solution illustrating that the diffusive term is 
well resolved. The test was carried out for various number of nodes n x = 32, 64, 128 and 288. Even with only 32 grid 
points the diffusion is resolved with an accuracy better than 1%. 

The tests illustrate the advantageous behaviour of the chosen combination of a WENO- scheme with a fourth-order 
discretization of the diffusive terms which resulted in a low numerical diffusion and small absolute errors for both 
modes of transport, pure convection and pure diffusion, respectively. 
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Figure 3: Comparison of numerical and analytical results of the diffusion concentration boundary layer 



4. Mesh sensitivity in 2D case 
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Figure 4: Schematic representation of the computational domain. The scatter shows every \O th grid line of the major grid used for 
the velocity field. The mesh for the scalar was further refined by factors ofR - 2 and R = 3. 

Fig. [4] shows a schematic representation of the computational domain. A quadratic 2D domain was chosen with 
an edge length of 5L. Periodic boundary conditions in the horizontal direction and a free slip condition at the top and 
bottom were applied. The base grid size was n x = 400 and n z = 256 in the x- and z-directions, respectively. The mesh 
was stretched in the z-direction to obtain a finer resolution near the top where a steep concentration gradient occurs. 
The procedure for the stretching is controlled by the parameter 6 = 3, the domain size is given by 



z(0) 
z(n z ) 




5 



(27) 
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Fork = l,...,n z - 1 



Subsequently we have, 



Z<p = (6/2)xk/n z (28) 
zi = 6/2 



co = (29) 
tanh(z\) 



The node distribution in the z-direction from z(0) = to z(n z ) = 5L is then given by, 

z(k) = (l-w)xz(0) + wxzW (30) 

for k = 1, ...,ft z - 1. 

The mesh sensitivity for resolving the flow and concentration fields were tested in several subgrid mesh refinement 
studies as described in the following. 



4.1. Flow-field 

To verify that the flow-field was fully resolved on the chosen 400 x 256 base mesh, the grid was refined in all 
directions by factors of 1.5 and 2, respectively. Fig. [5] shows the contour plots and the velocity profile obtained from 
the simulations with the base grid and the mesh refined by factor 2 (R = 2) (with 800 x 512 points) after t = 45 L/U. 
The contour plots of the flow-field using the refined mesh did not show any visible changes in the flow structures 
(Fig.[5b|). This is further confirmed by the vertical velocity profiles along a horizontal line at z = 4L . The profiles 
show a nice convergence verifying that the velocity field is fully resolved on the 400x256 grid which was subsequently 
used in all further cases. 
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4.2. Scalar transport: gas concentration field in buoyant flow 

As mentioned previously, we used a dual-mesh approach in which the scalar was resolved on a finer mesh than 
the one used for the velocity. Various levels of refinement were employed as illustrated in Fig [2] In this section, the 
mesh sensitivity for the scalar transport using this dual mesh approach is evaluated. 
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(a) Domain after t- 45 L/U on standard mesh 



(b) Domain after t- 45 L/U with refined submesh R = 3 




(c) Zoomed view after t - 45 L/U on standard mesh 



(d) Zoomed view after t- 45 L/U with refined submesh R = 3 



Figure 6: comparison of the gas concentration field after t - 45 L/U on standard mesh n x - 400 and n z - 256 and with a dual 
subgrid in place three times as fine (see Fig. [2c| . The gas concentration field is resolved in sharper detail with less smearing. 

Fig. [6] shows a comparison of the non-dimensional gas concentration contour plots that visualise the development 
of the scalar transport at t = 45 L/U using the base mesh (400 x 256) for both velocity and scalar and the dual mesh 
approach with refinement factor 3 applied to the scalar. The Schmidt number is S c = 500 which is equivalent to the 
diffusion of oxygen in water. The convective instability was triggered by adding random disturbances of 1% to the 
temperature field after letting it evolve for t = 10 L/U time units in order to avoid the triggering of the instability to 
depend on the mesh size or numerical round off error. The same disturbance field was used for all simulations. In 
general, both concentration fields in Figs. [6a| and [6b] reveal the same structures of downwards plumes. However, a 
zoomed view of the top region near the water surface reveals a more detailed representation of the gas concentration 
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field when the dual mesh approach is used (Fig. [6c]and[6d]). Please note that the grid plotted in those figures is the 
base grid and not the refined mesh used for the scalar transport. 





Figs. [7] and [8] show line plots of the scalar field at various locations within the domain. The locations are across 
or along the typical mushroom pattern that develops as a result of the convective instability where sharp gradients in 
the scalar field are present. Solving the scalar on the finer subgrid shows a significant improvement in resolution. The 
R - 2 refinement shows a big improvement in deeper regions where the base mesh is relatively coarse and the scalar 
distribution is maintained better. In the far field (z < 4) the scalar concentration profiles for the refined cases R = 2, 
R = 3 and R = 5 converge to nearly identical values (Fig.[8a|. 

The improved resolution becomes even more relevant when the spatially integrated total scalar concentration in 
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Figure 9: Comparison of the total non-dimensional scalar concentration and temperature T*over time for different levels of 
subgrid refinement 



the domain over time is considered. Fig. [9a| shows the total concentration over time for Sc = 500. Up to a time of 
t = 30 L/U the gas transfer is dominated by diffusion. Subsequently, the instability induces a convective flow that 
enhances the mass transfer. The typical mushroom patterns start penetrating the deeper regions of the domain. It is 
here where the refined submesh shows a much improved resolution with a continuous increase in the concentration 
levels whilst the standard mesh shows a drop in concentration levels. The drop occurs when the scalar reaches the 



region for z < 2.5L (where the mesh becomes significantly coarser) after around t = 55 L/U (Fig. [9a]). This points out 
an insufficient resolution of the scalar transport in this region. This effect was not present in the refined cases (Fig.[9a|. 



The same is found for the transport of the non-dimensionalized temperature T* (Fig. 9b ). The grid refinement study 
for the temperature transport shows a similar trend as seen for the concentration field in Fig. [9b] On the coarse mesh 
fluctuations become evident after t = 50 L/U whereas the refined cases do not exhibit such temperature fluctuations. 
Again the results are identical for both refined cases (R - 2, R - 3 and R - 5). 



5. Comparison to Experiments 

The oxygen concentration field was measured in experiments at KIT by a Laser Induced Fluorescence (LIF) 
method based on the oxygen quenching phenomenon developed by Vaughan & Weber EH . The method enables an 
instantaneous 2D visualisation with good resolution of the dissolved oxygen concentration l22l . The experiments 
were performed in a tank where the surface temperature was 3 °C lower than the bulk temperature of the water. The 
equivalent temperature boundary condition was applied in the numerical simulations. 



Fig. [10| shows a comparison of 2D-LIF images to our 2D DNS results. Note that the DNS results show the top 
section of the domain that has the same dimension as the LIF-maps. The actual experimental domain domain was 
much larger so the sides and bottom in these plots can be considered as open boundaries. For reasons of better 
comparison the timescale was set to t - L/U from the moment when the flow field started moving. Though the 
numerical results are only 2D, whereas the real problem is of course 3D, a very good qualitative agreement with 
the experiment is observed. Especially between t = 16 L/U and t = 28 L/U both the spatial distance between high 
concentration plumes and the size of the eddies were found to be similar in the experiment and the simulation. Because 
of the low diffusivity of oxygen in water and the rather low turbulent flow the plumes of high oxygen concentration 
retain their fine structures. This means that the steep concentration gradients do not smear out because of excessive 
numerical diffusion. As a result good qualitative agreement between the numerical simulations and the experimental 
data is obtained. 
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(b) Numerical Results 

Figure 10: Comparison of flow structures. High Oxygen concentration plumes of LIF measurements (Fig. \lOa\ t22i and DNS 
results (Fig.UUbi. In both cases the surface temperature is 3°C colder than the bulk temperature 



6. Conclusion 

To accurately resolve the mass transport for a scalar with low diffusivity on a stretched and staggered mesh, 
the classical WENO scheme of Liu et al iTTTTl was implemented to discretize the convective terms, while a fourth 
order accurate central discretization was used for the diffusion. The flow field was approximated with fourth-order 
accurate central discretizations for both convection and diffusion. Because the diffusivity of the scalars of interest 
is up to almost three orders of magnitude smaller than the molecular diffusivity of the ambient fluid, the resolution 
requirements for the transported scalar are much higher. Hence, to save computing time, a dual meshing approach 
was employed in which the scalar transport equations were discretized on a finer mesh than the flow field. The 
discretization of scalar convection and diffusion were tested in both ID and 2D cases. The ID tests showed that 
the spatial discretization of the convective term achieved a third-order accuracy on non-uniform meshes, while the 
discretization of the diffusive term was shown to achieve a very low relative error even on coarse meshes at a Schmidt 
number of Sc = 500. For the 2D case a combined active and passive scalar problem was simulated. The dual 
meshing approach showed a significant improvement in accuracy of the scalar field resolution even for a moderate 
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refinement by a factor of two. Subsequent refinements that were carried out using factors of up to five times, only 
showed marginal further improvements in accuracy. Additionally a qualitative comparison of the numerical results to 
experimentally visualized oxygen concentration fields in water showed similar structures below the air-water interface 
even though the numerical simulations were 2-dimensional. 
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